//#define BOOST_TEST_MODULE interp_twomesh tests
#include <boost/test/unit_test.hpp>
using boost::unit_test::test_suite;

#include <feel/options.hpp>
#include <feel/feelalg/backend.hpp>
#include <feel/feeldiscr/functionspace.hpp>
#include <feel/feelfilters/gmsh.hpp>
#include <feel/feelfilters/exporter.hpp>
#include <feel/feelvf/vf.hpp>
//#include <feel/feelfilters/gmshtensorizeddomain.hpp>
#include <feel/feelfilters/geotool.hpp>
#include <feel/feelcore/applicationxml.hpp>
#include <iostream>

using namespace Feel;
using namespace Feel::vf;



namespace test_interp_twomesh
{

    typedef Application Application_type;
    //typedef ApplicationXML Application_type;

    typedef boost::shared_ptr<Application_type> Application_ptrtype;

    enum type_champs
        {
            ScalarTest = 0,
            VectorialTest
        };


    /*_________________________________________________*
     * Options
     *_________________________________________________*/

    inline
    po::options_description
    makeOptions()
    {
        po::options_description desc_options("test_interp_twomesh options");
        desc_options.add_options()
            ("hsize1", po::value<double>()->default_value( 0.05 ), "mesh size1")
            ("hsize2", po::value<double>()->default_value( 0.03 ), "mesh size2")
            ("gmsh", po::value<float>()->default_value( 2.1 ), " version of gmsh(2.0 or 2.1)")
            ;
        return desc_options.add( Feel::feel_options() );
    }

    /*_________________________________________________*
     * About
     *_________________________________________________*/

    inline
    AboutData
    makeAbout()
    {
        AboutData about( "Test_Interp_Conv" ,
                         "Test_Interp_Conv" ,
                         "0.1",


                         "Test1 verify geomap such as phi*ph^{-1}=Id ; Test2: on gamma, ||..v(u1)-..v(u2)|| < epsilon",
                         Feel::AboutData::License_GPL,
                         "Copyright (c) 2010 Universite Joseph Fourier");

        about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
        return about;

    }

    /*_________________________________________________*
     * createMesh
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh( std::string __name,std::string __str)
    {
        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( new mesh_type );

        Gmsh gmsh;
        gmsh.setOrder(OrderGeo);

        std::string fname = gmsh.generate( __name, __str,false,false,false  );

        ImporterGmsh<mesh_type> import( fname );

        float gmsh_version = 2.1;//this->vm()["gmsh"].template as<float>();
        if (gmsh_version==2.0) import.setVersion( "2.0" );
        else if (gmsh_version==2.1) import.setVersion( "2.1" );

        mesh->accept( import );
        mesh->components().set ( MESH_RENUMBER|MESH_UPDATE_EDGES|MESH_UPDATE_FACES|MESH_CHECK );
        mesh->updateForUse();

        return mesh;

    }

    /*_________________________________________________*
     * createMesh_rect_p
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_rect_p( std::string __name,double meshSize)
    {
        GeoTool::Node x1(-1,-1);
        GeoTool::Node x2(1,1);
        GeoTool::Rectangle R1( meshSize,"LEFT",x1,x2);
        //R1.setMarker(_type="line",_name="Interface",_marker2=true);
        R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
        R1.setMarker(_type="surface",_name="Omega1",_marker1=true);


        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     R1.geoStr() ));

        return mesh;

    }

    /*_________________________________________________*
     * createMesh_rect_p1
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_rect_p1( std::string __name,double meshSize)
    {
        GeoTool::Node x1(-1,-1);
        GeoTool::Node x2(0.3,1);
        GeoTool::Rectangle R1( meshSize,"LEFT",x1,x2);
        R1.setMarker(_type="line",_name="Interface",_marker2=true);
        R1.setMarker(_type="line",_name="Bord",_marker1=true,_marker3=true,_marker4=true);
        R1.setMarker(_type="surface",_name="Omega1",_marker1=true);


        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     R1.geoStr() ));

        return mesh;

    }

    /*_________________________________________________*
     * createMesh_rect_p2
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_rect_p2( std::string __name,double meshSize)
    {
        GeoTool::Node x3(0.3,-1);
        GeoTool::Node x4(1,1);
        GeoTool::Rectangle R(meshSize,"RIGHT",x3,x4);
        R.setMarker(_type="line",_name="Interface",_marker4=true);
        R.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true);
        R.setMarker(_type="surface",_name="Omega2",_marker1=true);

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     R.geoStr() ));

        return mesh;
    }


    /*_________________________________________________*
     * createMesh
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_rectcirc_p1( std::string __name,double meshSize)
    {

        GeoTool::Node x21(-1,-1);
        GeoTool::Node x22(1,1);
        GeoTool::Rectangle R(meshSize,"RIGHT",x21,x22);
        R.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
        R.setMarker(_type="surface",_name="Omega3",_marker1=true);
	
        GeoTool::Node x_centre(0.0,0.0);
        GeoTool::Node x_bord(0.3,0.0);
        GeoTool::Circle C(meshSize/2,"DISQUE",x_centre,x_bord);
        C.setMarker(_type="line",_name="Interface",_marker1=true);
        C.setMarker(_type="surface",_name="Omega2",_marker1=true);

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     (R-C).geoStr() ));

        return mesh;

    }

    /*_________________________________________________*
     * createMesh
     *_________________________________________________*/

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_rectcirc_p2( std::string __name,double meshSize)
    {
        //GeoTool::Rectangle R4(meshSize2,"RIGHT",x21,x22);
        //R4.setMarker(_type="line",_name="Bord",_marker1=true,_marker2=true,_marker3=true,_marker4=true);
        //R4.setMarker(_type="surface",_name="Omega3",_marker1=true);
        GeoTool::Node x_centre(0.0,0.0);
        GeoTool::Node x_bord(0.3,0.0);
	
        GeoTool::Circle C(meshSize,"DISQUE",x_centre,x_bord);
        C.setMarker(_type="line",_name="Interface",_marker1=true);
        C.setMarker(_type="surface",_name="Omega2",_marker1=true);

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     C.geoStr() ));

        return mesh;


    }

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_spec_p1( std::string __name,double meshSize)
    {


        GeoTool::Node xbid(-1,-1);
        GeoTool::Special_1a S1a(meshSize,"spec",xbid);
        S1a.setMarker(_type="line",_name="inner",_marker1=true);
        S1a.setMarker(_type="line",_name="outer",_marker2=true);
        S1a.setMarker(_type="line",_name="boundLeft",_marker3=true);
        S1a.setMarker(_type="line",_name="boundRight",_marker4=true);
        S1a.setMarker(_type="surface",_name="membrane",_marker1=true);

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     S1a.geoStr() ));

        return mesh;

    }

    template<uint Dim,uint OrderGeo>
    boost::shared_ptr< Mesh<Simplex<Dim,OrderGeo,Dim> > >
    createMesh_spec_p2( std::string __name,double meshSize)
    {


        GeoTool::Node xbid(-1,-1);
        GeoTool::Special_1b S1b(meshSize,"spec2",xbid);
        S1b.setMarker(_type="line",_name="interface",_marker1=true);
        S1b.setMarker(_type="line",_name="in",_marker2=true);
        S1b.setMarker(_type="line",_name="out",_marker3=true);
        S1b.setMarker(_type="surface",_name="fluide",_marker1=true);

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<mesh_type> mesh_ptrtype;
        mesh_ptrtype mesh( createMesh<Dim,OrderGeo>( __name + mesh_type::shape_type::name(),
                                                     S1b.geoStr() ));

        return mesh;

    }




    template<typename element_type>
    void
    projectionH1(element_type & u,Application_ptrtype & test_app)
    {
        typedef Backend<double> backend_type;
        typedef boost::shared_ptr<backend_type> backend_ptrtype;
        typedef backend_type::sparse_matrix_type sparse_matrix_type;
        typedef backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
        typedef backend_type::vector_type vector_type;
        typedef backend_type::vector_ptrtype vector_ptrtype;

        typedef typename element_type::functionspace_type functionspace_type;
        typedef typename element_type::functionspace_ptrtype functionspace_ptrtype;

        static const uint16_type OrderChamp = functionspace_type::basis_type::nOrder;
        static const uint16_type OrderGeo = 2;

        functionspace_ptrtype Xh=u.functionSpace();

        backend_ptrtype backend(backend_type::build( test_app->vm() ));

        sparse_matrix_ptrtype A = backend->newMatrix(Xh,Xh);
        vector_ptrtype F = backend->newVector(Xh);

        AUTO( pi, M_PI );
        AUTO( Uex, cos(pi*Px()/2.)*cos(pi*Py()/2.) );
        AUTO( grad_Uex , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );


        form2( Xh, Xh, A, _init=true ) =
            integrate( elements(Xh->mesh()),_Q<2*(OrderChamp*OrderGeo)>(),
                                              idt(u)*id(u)
                                              +gradt(u)*trans(grad(u)) );
        form1( Xh, F, _init=true ) = integrate( elements(Xh->mesh()),
                                                _Q<18>(),
                                                Uex*id(u)+grad(u)*grad_Uex );

                       //form1( Xh, F) += integrate( markedfaces(mesh, mesh->markerName("Interface")),trans(grad_Uex)*N()*id(u) );
        A->close();
        F->close();

        /*
          form2( Xh, Xh, A ) +=
          on( boundaryfaces(Xh->mesh()) , u, F, Uex );
        */

        //      on( markedfaces(mesh1, mesh1->markerName("Bord")) , u, F, cst(0.) );

        vector_ptrtype X= backend->newVector( Xh );
        backend->solve( A, A, X, F, false );
        u = *X;
    }

        /*_________________________________________________*
         * run_test_element
         *_________________________________________________*/



        template<uint Dim,uint OrderChamp,uint OrderGeo>
        void
        run_test_element(Application_ptrtype & test_app)
        {

            Log() << "[testGeoMap] start test geomap inverse\n";

            //-----------------------------------------------------------------------------------//

            typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
            typedef boost::shared_ptr<  mesh_type > mesh_ptrtype;

            typedef bases<Lagrange<OrderChamp,Scalar,PointSetFekete> > basis_type;
            typedef FunctionSpace<mesh_type, basis_type> space_type;
            typedef boost::shared_ptr<space_type> space_ptrtype;
            typedef typename space_type::element_type element_type;

            //-----------------------------------------------------------------------------------//

            double meshSize1;
            double meshSize2;

            std::vector<boost::tuple<double,double> > meshSizeVec(5);
            //meshSizeVec[0]=boost::make_tuple(0.1,2);
            meshSizeVec[0]=boost::make_tuple(0.1/2.,2/2.);
            meshSizeVec[1]=boost::make_tuple(0.1/4.,2/4.);
            meshSizeVec[2]=boost::make_tuple(0.1/2.,2/8.);
            meshSizeVec[3]=boost::make_tuple(0.1/2.,2/16.);
            meshSizeVec[4]=boost::make_tuple(0.1/2.,2/32.);
            //meshSizeVec[5]=boost::make_tuple(0.1/2.,0.07/32.);

    
            std::ostringstream __ostr;
            __ostr<<"graph_elt_conv_"<<OrderGeo<<".data";


            std::ofstream out;
            //out.open("graph_conv.data"+__ostr.str(), std::ios::out);
            out.open(__ostr.str().c_str(), std::ios::out);

    
            std::vector<boost::tuple<double,double> >::const_iterator itmeshSize = meshSizeVec.begin();
            std::vector<boost::tuple<double,double> >::const_iterator itmeshSize_end = meshSizeVec.end();
            for ( ; itmeshSize!=itmeshSize_end; ++itmeshSize)
                {

                    meshSize1=2/32.;//32::get<0>(*itmeshSize);
                    meshSize2=boost::get<1>(*itmeshSize);

                    //-----------------------------------------------------------------------------------//
                    //Geometry


                    //-----------------------------------------------------------------------------------//
	

                    //-----------------------------------------------------------------------------------//

                    //Mesh

                    //std::ostringstream __ostr;__ostr<<
                    std::ostringstream __ostr_name1;
                    __ostr_name1 << "mesh_rect1_h" << meshSize1 <<"_";
                    std::ostringstream __ostr_name2;
                    __ostr_name2 << "mesh_rect2_h" << meshSize2 <<"_";
                    std::ostringstream __ostr_name3;
                    __ostr_name3 << "mesh_rectcirc1_h" << meshSize1 <<"_";
                    std::ostringstream __ostr_name4;
                    __ostr_name4 << "mesh_rectcirc2_h" << meshSize2 <<"_";


                    mesh_ptrtype mesh1 = createMesh_rect_p<Dim,OrderGeo>( __ostr_name1.str()/* + mesh_type::shape_type::name()*/,
                                                                          meshSize1 );
                    mesh_ptrtype mesh2 = createMesh_rect_p<Dim,OrderGeo>( __ostr_name2.str()/* + mesh_type::shape_type::name()*/,
                                                                          meshSize2 );
                    mesh_ptrtype mesh3 = createMesh_rectcirc_p1<Dim,OrderGeo>(__ostr_name3.str() /*+ mesh_type::shape_type::name()*/,
                                                                              meshSize1 );
                    mesh_ptrtype mesh4 = createMesh_rectcirc_p1<Dim,OrderGeo>( __ostr_name4.str() /*+ mesh_type::shape_type::name()*/,
                                                                               meshSize2 );

                    space_ptrtype Xh1 = space_type::New( mesh1 );
                    space_ptrtype Xh2 = space_type::New( mesh2 );
                    space_ptrtype Xh3 = space_type::New( mesh3 );
                    space_ptrtype Xh4 = space_type::New( mesh4 );

                    element_type u1( Xh1, "u1" );
                    element_type u2( Xh2, "u2" );
                    element_type u3( Xh3, "u3" );
                    element_type u4( Xh4, "u4" );

                    //-----------------------------------------------------------------------------------//

                    //AUTO (e ,exp(Px()*Py())*sin(2*M_PI*Px()));
                    //AUTO (f , sin(2*M_PI*Px())*sin(2*M_PI*Py()));
                    AUTO (h , cos(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
                    //AUTO (g , Px()*(Px()-1)*Py()*(Py()-1) );
	
                    //AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );

                    u1 = vf::project( Xh1, elements(mesh1), h );
                    u2 = vf::project( Xh2, elements(mesh2), h );
                    u3 = vf::project( Xh3, elements(mesh3), h );
                    u4 = vf::project( Xh4, elements(mesh4), h );

                    //-----------------------------------------------------------------------------------//

                    double  __errId1 = std::sqrt(integrate( elements(mesh1),_Q<2*(OrderChamp+OrderGeo-1)>(),
                                                            (idv(u2)-h/*idv(u4)*/)*(idv(u2)-h/*idv(u4)*/) ).evaluate()(0,0));
                    double  __errId2 = std::sqrt(integrate( elements(mesh3),_Q<2*(OrderChamp+OrderGeo-1)>(),
                                                            (idv(u4)-h/*idv(u4)*/)*(idv(u4)-h/*idv(u4)*/) ).evaluate()(0,0));

                    std::cout << "Error L2 : " << meshSize2 << " " << __errId1 << " " << __errId2 << "\n";
                    out << meshSize2 << " " << __errId1 << " " << __errId2 << "\n";
                }
            //      }
            Log() << "[testGeoMap] finish test geomap inverse\n";

        } //end run


	//-----------------------------------------------------------------------------------//
	//-----------------------------------------------------------------------------------//
	//-----------------------------------------------------------------------------------//
	//-----------------------------------------------------------------------------------//
	//-----------------------------------------------------------------------------------//


    template<uint Dim,uint OrderChamp,uint OrderGeo>
    void
    run_test_convergence_h(Application_ptrtype & test_app)
    {

        Log() << "[testGeoMap] start test geomap inverse\n";

        //-----------------------------------------------------------------------------------//

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<  mesh_type > mesh_ptrtype;

        typedef bases<Lagrange<OrderChamp,Scalar,PointSetFekete> > basis_type;
        typedef FunctionSpace<mesh_type, basis_type> space_type;
        typedef boost::shared_ptr<space_type> space_ptrtype;
        typedef typename space_type::element_type element_type;

        //-----------------------------------------------------------------------------------//

        double meshSize1;
        double meshSize2;


        std::vector<boost::tuple<double,double> > meshSizeVec(8);
        meshSizeVec[0]=boost::make_tuple(0.1,2.);
        meshSizeVec[1]=boost::make_tuple(0.1,1.);
        meshSizeVec[2]=boost::make_tuple(0.1,1/2.);
        meshSizeVec[3]=boost::make_tuple(0.1/2.,1/4.);
        meshSizeVec[4]=boost::make_tuple(0.1/4.,1/8.);
        meshSizeVec[5]=boost::make_tuple(0.1/4.,1/12.);
        meshSizeVec[6]=boost::make_tuple(0.1/2.,1/16.);
        meshSizeVec[7]=boost::make_tuple(0.1/2.,1/24.);
        //meshSizeVec[6]=boost::make_tuple(0.1/2.,1/32.);
        //meshSizeVec[7]=boost::make_tuple(0.1/2.,1/48.);
        //meshSizeVec[8]=boost::make_tuple(0.1/2.,1/64.);

        //std::ostringstream __ostr;
        //__ostr<<"graph_elt_conv_"<<OrderGeo<<".data";

        std::ostringstream __ostrrr;
        __ostrrr<<"graph_boundary_conv_h_"<<OrderGeo<<".data";
        std::ofstream out;
        out.open(__ostrrr.str().c_str(), std::ios::out);

        std::vector<boost::tuple<double,double> >::const_iterator itmeshSize = meshSizeVec.begin();
        std::vector<boost::tuple<double,double> >::const_iterator itmeshSize_end = meshSizeVec.end();
        for ( ; itmeshSize!=itmeshSize_end; ++itmeshSize)
            {

                meshSize2=1./24.;//32.;//256.;//32::get<0>(*itmeshSize);
                meshSize1=boost::get<1>(*itmeshSize);

                //-----------------------------------------------------------------------------------//

                //Mesh

                std::string __nameTest="boundary";
                //std::ostringstream __ostr;__ostr<<
                std::ostringstream __ostr_name1;
                __ostr_name1 << __nameTest << "_mesh_rect1_h" << meshSize1 <<"_";
                std::ostringstream __ostr_name2;
                __ostr_name2 << __nameTest << "mesh_rect2_h" << meshSize2 <<"_";
                std::ostringstream __ostr_name3;
                __ostr_name3 << __nameTest << "mesh_rectcirc1_h" << meshSize1 <<"_";
                std::ostringstream __ostr_name4;
                __ostr_name4 << __nameTest << "mesh_rectcirc2_h" << meshSize2 <<"_";

                mesh_ptrtype mesh1 = createMesh_rect_p1<Dim,OrderGeo>( __ostr_name1.str(), meshSize1 );
                mesh_ptrtype mesh2 = createMesh_rect_p2<Dim,OrderGeo>( __ostr_name2.str(), meshSize2 );
                //mesh_ptrtype mesh1 = createMesh_spec_p1<Dim,OrderGeo>( __ostr_name1.str(), meshSize1 );
                //mesh_ptrtype mesh2 = createMesh_spec_p2<Dim,OrderGeo>( __ostr_name2.str(), meshSize2 );

                mesh_ptrtype mesh3 = createMesh_rectcirc_p1<Dim,OrderGeo>(__ostr_name3.str(), meshSize1 );
                mesh_ptrtype mesh4 = createMesh_rectcirc_p2<Dim,OrderGeo>( __ostr_name4.str(), meshSize2/6. );

                //-----------------------------------------------------------------------------------//

                space_ptrtype Xh1 = space_type::New( mesh1 );uint nb1dof=Xh1->nDof();
                space_ptrtype Xh2 = space_type::New( mesh2 );uint nb2dof=Xh2->nDof();
                space_ptrtype Xh3 = space_type::New( mesh3 );uint nb3dof=Xh3->nDof();
                space_ptrtype Xh4 = space_type::New( mesh4 );uint nb4dof=Xh4->nDof();

                element_type u1( Xh1, "u1" );
                element_type u2( Xh2, "u2" );
                element_type u3( Xh3, "u3" );
                element_type u4( Xh4, "u4" );

                //-----------------------------------------------------------------------------------//
                /*
                  AUTO (h , cos(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
                  AUTO (dxh , -0.5*M_PI*sin(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
                  AUTO (dyh , -0.5*M_PI*cos(0.5*M_PI*Px())*(sin(0.5*M_PI*Py())));
                  AUTO (grad_h, vec(dxh,dyh));
                */
                AUTO( pi, M_PI );
                AUTO( h, cos(pi*Px()/2.)*cos(pi*Py()/2.) );
                AUTO( grad_h , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );

                u1 = vf::project( Xh1, elements(mesh1), h );
                u2 = vf::project( Xh2, elements(mesh2), h );

                projectionH1<element_type>(u3,test_app);
                projectionH1<element_type>(u4,test_app);




                //-----------------------------------------------------------------------------------//

                double  __errId1 = integrate( markedfaces(mesh2, mesh2->markerName("Interface")),
                                              _Q<2*(OrderChamp*OrderGeo)>(),//_Q<2*(OrderChamp+OrderGeo-1)>(),
                                              (idv(u1)-h)*(idv(u1)-h) ).evaluate()(0,0);

                double  __errId2 = integrate( markedfaces(mesh4, mesh4->markerName("Interface")),
                                              _Q<2*(OrderChamp*OrderGeo)>(), //_Q<2*(OrderChamp+OrderGeo-1)>(),
                                              (idv(u3)-h)*(idv(u3)-h) ).evaluate()(0,0);

                double  __errGrad1 = integrate( markedfaces(mesh2, mesh2->markerName("Interface")),
                                                _Q<2*(OrderChamp*OrderGeo)>(),//_Q<2*(OrderChamp+OrderGeo-1)>(),
                                                (gradv(u1)-trans(grad_h))*(trans(gradv(u1))-grad_h) ).evaluate()(0,0);

                double  __errGrad2 = integrate( markedfaces(mesh4, mesh4->markerName("Interface")),
                                                _Q<2*(OrderChamp*OrderGeo)>(),//_Q<2*(OrderChamp+OrderGeo-1)>(),
                                                (gradv(u3)-trans(grad_h))*(trans(gradv(u3))-grad_h) ).evaluate()(0,0);

                //-----------------------------------------------------------------------------------//

                std::cout << "Error L2 : " << meshSize2 << " " << __errId2 << "\n";
                std::cout << "Error H1 : " << meshSize2 << " " << __errId2+__errGrad2 << "\n";

                //-----------------------------------------------------------------------------------//

                out << nb1dof << " "
                     << std::sqrt(__errId1) << " "
                     << std::sqrt(__errId1+__errGrad1) << " "
                     << nb3dof << " "
                     << std::sqrt(__errId2) << " "
                     << std::sqrt(__errId2+__errGrad2) << "\n";

            }


    } //end run





    template<uint Dim,uint OrderChamp>
    void
    run_test_boundary2(Application_ptrtype & test_app)
    {

        Log() << "[testGeoMap] start test boundary2\n";

        //-----------------------------------------------------------------------------------//

        const uint OrderGeo = 5;

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_base_type;
        typedef Mesh<Simplex<Dim,1,Dim> > mesh_1_type;
        typedef Mesh<Simplex<Dim,2,Dim> > mesh_2_type;
        typedef Mesh<Simplex<Dim,3,Dim> > mesh_3_type;
        typedef Mesh<Simplex<Dim,4,Dim> > mesh_4_type;
        typedef Mesh<Simplex<Dim,5,Dim> > mesh_5_type;

        typedef boost::shared_ptr< mesh_base_type > mesh_base_ptrtype;
        typedef boost::shared_ptr< mesh_1_type > mesh_1_ptrtype;
        typedef boost::shared_ptr< mesh_2_type > mesh_2_ptrtype;
        typedef boost::shared_ptr< mesh_3_type > mesh_3_ptrtype;
        typedef boost::shared_ptr< mesh_4_type > mesh_4_ptrtype;
        typedef boost::shared_ptr< mesh_5_type > mesh_5_ptrtype;

        typedef bases<Lagrange<6,Scalar,PointSetFekete> > basis_1_type;
        typedef bases<Lagrange<6,Scalar,PointSetFekete> > basis_2_type;
        typedef bases<Lagrange<6,Scalar,PointSetFekete> > basis_3_type;
        typedef bases<Lagrange<6,Scalar,PointSetFekete> > basis_4_type;
        typedef bases<Lagrange<6,Scalar,PointSetFekete> > basis_5_type;

        typedef FunctionSpace<mesh_base_type, basis_5_type> space_base_type;
        typedef FunctionSpace<mesh_1_type, basis_1_type> space_1_type;
        typedef FunctionSpace<mesh_2_type, basis_2_type> space_2_type;
        typedef FunctionSpace<mesh_3_type, basis_3_type> space_3_type;
        typedef FunctionSpace<mesh_4_type, basis_4_type> space_4_type;
        typedef FunctionSpace<mesh_5_type, basis_5_type> space_5_type;

        typedef boost::shared_ptr<space_base_type> space_base_ptrtype;
        typedef boost::shared_ptr<space_1_type> space_1_ptrtype;
        typedef boost::shared_ptr<space_2_type> space_2_ptrtype;
        typedef boost::shared_ptr<space_3_type> space_3_ptrtype;
        typedef boost::shared_ptr<space_4_type> space_4_ptrtype;
        typedef boost::shared_ptr<space_5_type> space_5_ptrtype;

        typedef typename space_base_type::element_type element_base_type;
        typedef typename space_1_type::element_type element_1_type;
        typedef typename space_2_type::element_type element_2_type;
        typedef typename space_3_type::element_type element_3_type;
        typedef typename space_4_type::element_type element_4_type;
        typedef typename space_5_type::element_type element_5_type;

        //-----------------------------------------------------------------------------------//

        double meshSize1;
        double meshSize2;

        std::vector<boost::tuple<double,double> > meshSizeVec(5);
        meshSizeVec[0]=boost::make_tuple(0.1,1/2.);
        meshSizeVec[1]=boost::make_tuple(0.1/2.,1/4.);
        meshSizeVec[2]=boost::make_tuple(0.1/4.,1/8.);
        meshSizeVec[3]=boost::make_tuple(0.1/2.,1/16.);
        meshSizeVec[4]=boost::make_tuple(0.1/2.,1/32.);
        //meshSizeVec[5]=boost::make_tuple(0.1/2.,0.07/32.);

    
        std::ostringstream __ostr;
        //__ostr<<"graph_elt_convTOT_"<<OrderGeo<<".data";
        __ostr<<"graph_elt_convTOT.data";

        std::ofstream out;
        //out.open("graph_conv.data"+__ostr.str(), std::ios::out);
        out.open(__ostr.str().c_str(), std::ios::out);

    
        std::vector<boost::tuple<double,double> >::const_iterator itmeshSize = meshSizeVec.begin();
        std::vector<boost::tuple<double,double> >::const_iterator itmeshSize_end = meshSizeVec.end();
        for ( ; itmeshSize!=itmeshSize_end; ++itmeshSize)
            {

                meshSize1=2/32.;//32::get<0>(*itmeshSize);
                meshSize2=boost::get<1>(*itmeshSize);

                //-----------------------------------------------------------------------------------//

                //Mesh

                std::string __nameTest="boundary2";
                //std::ostringstream __ostr;__ostr<<

                std::ostringstream __ostr_base_name1;
                __ostr_base_name1 << __nameTest << "_mesh_base_rect1_h" << meshSize1 <<"_";
                std::ostringstream __ostr_base_name3;
                __ostr_base_name3 << __nameTest << "_mesh_base_rectcirc1_h" << meshSize1 <<"_";

                std::ostringstream __ostr_name1;
                __ostr_name1 << __nameTest << "_mesh_rect1_h" << meshSize1 <<"_";
                std::ostringstream __ostr_name2;
                __ostr_name2 << __nameTest << "mesh_rect2_h" << meshSize2 <<"_";
                std::ostringstream __ostr_name3;
                __ostr_name3 << __nameTest << "mesh_rectcirc1_h" << meshSize1 <<"_";
                std::ostringstream __ostr_name4;
                __ostr_name4 << __nameTest << "mesh_rectcirc2_h" << meshSize2 <<"_";
	
                //mesh_base_ptrtype mesh1 = createMesh_rect_p1<Dim,mesh_base_type::nOrder>( __ostr_base_name1.str(), meshSize1 );
                mesh_base_ptrtype mesh1 = createMesh_spec_p1<Dim,mesh_base_type::nOrder>( __ostr_base_name1.str(), meshSize1 );

                mesh_base_ptrtype mesh3 = createMesh_rectcirc_p1<Dim,mesh_base_type::nOrder>( __ostr_base_name3.str(), meshSize1 );
                /*
                  mesh_1_ptrtype mesh2_1 = createMesh_rect_p2<Dim,1>( __ostr_name2.str()+"_1", meshSize2 );
                  mesh_2_ptrtype mesh2_2 = createMesh_rect_p2<Dim,2>( __ostr_name2.str()+"_2", meshSize2 );
                  mesh_3_ptrtype mesh2_3 = createMesh_rect_p2<Dim,3>( __ostr_name2.str()+"_3", meshSize2 );
                  mesh_4_ptrtype mesh2_4 = createMesh_rect_p2<Dim,4>( __ostr_name2.str()+"_4", meshSize2 );
                  mesh_5_ptrtype mesh2_5 = createMesh_rect_p2<Dim,5>( __ostr_name2.str()+"_5", meshSize2 );
                */
                mesh_1_ptrtype mesh2_1 = createMesh_spec_p2<Dim,1>( __ostr_name2.str()+"_1", meshSize2 );
                mesh_2_ptrtype mesh2_2 = createMesh_spec_p2<Dim,2>( __ostr_name2.str()+"_2", meshSize2 );
                mesh_3_ptrtype mesh2_3 = createMesh_spec_p2<Dim,3>( __ostr_name2.str()+"_3", meshSize2 );
                mesh_4_ptrtype mesh2_4 = createMesh_spec_p2<Dim,4>( __ostr_name2.str()+"_4", meshSize2 );
                mesh_5_ptrtype mesh2_5 = createMesh_spec_p2<Dim,5>( __ostr_name2.str()+"_5", meshSize2 );


                mesh_1_ptrtype mesh4_1 = createMesh_rectcirc_p2<Dim,1>( __ostr_name4.str()+"_1", meshSize2/6.);
                mesh_2_ptrtype mesh4_2 = createMesh_rectcirc_p2<Dim,2>( __ostr_name4.str()+"_2", meshSize2/6.);
                mesh_3_ptrtype mesh4_3 = createMesh_rectcirc_p2<Dim,3>( __ostr_name4.str()+"_3", meshSize2/6.);
                mesh_4_ptrtype mesh4_4 = createMesh_rectcirc_p2<Dim,4>( __ostr_name4.str()+"_4", meshSize2/6.);
                mesh_5_ptrtype mesh4_5 = createMesh_rectcirc_p2<Dim,5>( __ostr_name4.str()+"_5", meshSize2/6.);
	
                //-----------------------------------------------------------------------------------//

                space_base_ptrtype Xh1 = space_base_type::New( mesh1 );
                space_base_ptrtype Xh3 = space_base_type::New( mesh3 );

                space_1_ptrtype Xh2_1 = space_1_type::New( mesh2_1 );
                space_2_ptrtype Xh2_2 = space_2_type::New( mesh2_2 );
                space_3_ptrtype Xh2_3 = space_3_type::New( mesh2_3 );
                space_4_ptrtype Xh2_4 = space_4_type::New( mesh2_4 );
                space_5_ptrtype Xh2_5 = space_5_type::New( mesh2_5 );

                space_1_ptrtype Xh4_1 = space_1_type::New( mesh4_1 );
                space_2_ptrtype Xh4_2 = space_2_type::New( mesh4_2 );
                space_3_ptrtype Xh4_3 = space_3_type::New( mesh4_3 );
                space_4_ptrtype Xh4_4 = space_4_type::New( mesh4_4 );
                space_5_ptrtype Xh4_5 = space_5_type::New( mesh4_5 );

                element_base_type u1( Xh1, "u1" );
                element_base_type u3( Xh3, "u3" );

                element_1_type u2_1( Xh2_1, "u2_1" );
                element_2_type u2_2( Xh2_2, "u2_2" );
                element_3_type u2_3( Xh2_3, "u2_3" );
                element_4_type u2_4( Xh2_4, "u2_4" );
                element_5_type u2_5( Xh2_5, "u2_5" );

                element_1_type u4_1( Xh4_1, "u4_1" );
                element_2_type u4_2( Xh4_2, "u4_2" );
                element_3_type u4_3( Xh4_3, "u4_3" );
                element_4_type u4_4( Xh4_4, "u4_4" );
                element_5_type u4_5( Xh4_5, "u4_5" );

	
                //-----------------------------------------------------------------------------------//

                //AUTO (e ,exp(Px()*Py())*sin(2*M_PI*Px()));
                //AUTO (f , sin(2*M_PI*Px())*sin(2*M_PI*Py()));
                //AUTO (h , cos(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));


                AUTO (dxh , -0.5*M_PI*sin(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
                AUTO (dyh , -0.5*M_PI*cos(0.5*M_PI*Px())*(sin(0.5*M_PI*Py())));
                AUTO (grad_h, vec(dxh,dyh));
                //AUTO (g , Px()*(Px()-1)*Py()*(Py()-1) );
	
                AUTO (h , (Px()+1)*(Px()-1)*(Py()+1)*(Py()-1) );

                u1 = vf::project( Xh1, elements(mesh1), h );
                u3 = vf::project( Xh3, elements(mesh3), h );


                u2_1 = vf::project( Xh2_1, elements(mesh2_1), h );
                u2_2 = vf::project( Xh2_2, elements(mesh2_2), h );
                u2_3 = vf::project( Xh2_3, elements(mesh2_3), h );
                u2_4 = vf::project( Xh2_4, elements(mesh2_4), h );
                u2_5 = vf::project( Xh2_5, elements(mesh2_5), h );

                u4_1 = vf::project( Xh4_1, elements(mesh4_1), h );
                u4_2 = vf::project( Xh4_2, elements(mesh4_2), h );
                u4_3 = vf::project( Xh4_3, elements(mesh4_3), h );
                u4_4 = vf::project( Xh4_4, elements(mesh4_4), h );
                u4_5 = vf::project( Xh4_5, elements(mesh4_5), h );

	
                //-----------------------------------------------------------------------------------//
                /*
                  double  __errId1_1 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
                  _Q<2*(OrderChamp+1-1)>(),
                  (idv(u2_1)-h)*(idv(u2_1)-h) ).evaluate()(0,0);
                  double  __errId1_2 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
                  _Q<2*(OrderChamp+2-1)>(),
                  (idv(u2_2)-h)*(idv(u2_2)-h) ).evaluate()(0,0);
                  double  __errId1_3 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
                  _Q<2*(OrderChamp+3-1)>(),
                  (idv(u2_3)-h)*(idv(u2_3)-h) ).evaluate()(0,0);
                  double  __errId1_4 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
                  _Q<2*(OrderChamp+4-1)>(),
                  (idv(u2_4)-h)*(idv(u2_4)-h) ).evaluate()(0,0);
                  double  __errId1_5 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
                  _Q<2*(OrderChamp+5-1)>(),
                  (idv(u2_5)-h)*(idv(u2_5)-h) ).evaluate()(0,0);
                */
                double  __errId1_1 = integrate( markedfaces(mesh1, mesh1->markerName("inner")),
                                                _Q<2*(OrderChamp+5-1)>(),
                                                (idv(u2_1)-h)*(idv(u2_1)-h) ).evaluate()(0,0);
                double  __errId1_2 = integrate( markedfaces(mesh1, mesh1->markerName("inner")),
                                                _Q<2*(OrderChamp+5-1)>(),
                                                (idv(u2_2)-h)*(idv(u2_2)-h) ).evaluate()(0,0);
                double  __errId1_3 = integrate( markedfaces(mesh1, mesh1->markerName("inner")),
                                                _Q<2*(OrderChamp+5-1)>(),
                                                (idv(u2_3)-h)*(idv(u2_3)-h) ).evaluate()(0,0);
                double  __errId1_4 = integrate( markedfaces(mesh1, mesh1->markerName("inner")),
                                                _Q<2*(OrderChamp+5-1)>(),
                                                (idv(u2_4)-h)*(idv(u2_4)-h) ).evaluate()(0,0);
                double  __errId1_5 = integrate( markedfaces(mesh1, mesh1->markerName("inner")),
                                                _Q<2*(OrderChamp+5-1)>(),
                                                (idv(u2_5)-h)*(idv(u2_5)-h) ).evaluate()(0,0);


                double  __errId2_1 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                                                _Q<2*(OrderChamp+1-1)>(),
                                                (idv(u4_1)-h)*(idv(u4_1)-h) ).evaluate()(0,0);
                double  __errId2_2 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                                                _Q<2*(OrderChamp+2-1)>(),
                                                (idv(u4_2)-h)*(idv(u4_2)-h) ).evaluate()(0,0);
                double  __errId2_3 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                                                _Q<2*(OrderChamp+3-1)>(),
                                                (idv(u4_3)-h)*(idv(u4_3)-h) ).evaluate()(0,0);
                double  __errId2_4 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                                                _Q<2*(OrderChamp+4-1)>(),
                                                (idv(u4_4)-h)*(idv(u4_4)-h) ).evaluate()(0,0);
                double  __errId2_5 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                                                _Q<2*(OrderChamp+5-1)>(),
                                                (idv(u4_5)-h)*(idv(u4_5)-h) ).evaluate()(0,0);


                /*
                  double  __errId2 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                  _Q<2*(OrderChamp+OrderGeo-1)>(),
                  (idv(u4)-h)*(idv(u4)-h) ).evaluate()(0,0);

                  double  __errGrad1 = integrate( markedfaces(mesh1, mesh1->markerName("Interface")),
                  _Q<2*(OrderChamp+OrderGeo-1)>(),
                  (gradv(u2)-trans(grad_h))*(trans(gradv(u2))-grad_h) ).evaluate()(0,0);
                  double  __errGrad2 = integrate( markedfaces(mesh3, mesh3->markerName("Interface")),
                  _Q<2*(OrderChamp+OrderGeo-1)>(),
                  (gradv(u4)-trans(grad_h))*(trans(gradv(u4))-grad_h) ).evaluate()(0,0);


                  std::cout << "Error L2 : " << meshSize2 << " " << __errId1 << " " << meshSize2/6. <<" "<< __errId2 << "\n";
                  std::cout << "Error H1 : " << meshSize2 << " " << __errId1+__errGrad1 << " " << meshSize2/6. <<" "<< __errId2+__errGrad2 << "\n";
                */
                out << meshSize2 << " " 
                    << std::sqrt(__errId1_1) << " "
                    << std::sqrt(__errId1_2) << " "
                    << std::sqrt(__errId1_3) << " "
                    << std::sqrt(__errId1_4) << " "
                    << std::sqrt(__errId1_5) << " "
                    //<< std::sqrt(__errId1+__errGrad1) << " "
                    << meshSize2/6. << " "
                    << std::sqrt(__errId2_1) << " " 
                    << std::sqrt(__errId2_2) << " " 
                    << std::sqrt(__errId2_3) << " " 
                    << std::sqrt(__errId2_4) << " " 
                    << std::sqrt(__errId2_5) << " "
                    //<< std::sqrt(__errId2+__errGrad2) 
                    << "\n";
            }

        Log() << "[test_boundary] finish convergence test on boundary \n";

    } //end run

    template<typename element_type,typename mesh_ptrtype,uint OrderGeo>
    void
    integration(element_type &u,mesh_ptrtype & mesh, std::ofstream &out)
    {    
        const uint16_type nOrder = element_type::functionspace_type::basis_type::nOrder;

        AUTO( pi, M_PI );
        AUTO( h, cos(pi*Px()/2.)*cos(pi*Py()/2.) );
        AUTO( grad_h , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );

        double  __errId = integrate( markedfaces(mesh, mesh->markerName("Interface")),
                                     _Q<2*(nOrder+OrderGeo-1)>(),
                                     (idv(u)-h)*(idv(u)-h) ).evaluate()(0,0);
        double  __errGrad = integrate( markedfaces(mesh, mesh->markerName("Interface")),
                                       _Q<2*(nOrder+OrderGeo-1)>(),
                                       (gradv(u)-trans(grad_h))*(trans(gradv(u))-grad_h) ).evaluate()(0,0);

        out << u.functionSpace()->nDof() << " "  << std::sqrt(__errId) << " " << std::sqrt(__errId+__errGrad) << "\n";

    }
    template<uint Dim,uint OrderChamp,uint OrderGeo>
    void
    run_test_convergence_p(Application_ptrtype & test_app)
    {

        //if (test_app->preProcessing() == RUN_EXIT) return;
        //test_app->preProcessing();

        //Parameter p;
        //p=Parameter(_name="p",_type=CONT_ATTR,_cmdName="polynomial order",_values="1:1:1" );
        //p=Parameter(_name="p",_type=CONT_ATTR,_cmdName="p",_values="1:1:1" );
        //test_app->addParameter( Parameter(_name="dim",_type=DISC_ATTR,_values=boost::lexical_cast<std::string>( Dim  ).c_str()) );

        //test_app->addParameter( Parameter(_name="dim",_type=DISC_ATTR,_values="2") );

        //test_app->addParameter(p);
        /*
          std::vector<Parameter> depend;
          std::vector<std::string> funcs;
          depend.push_back(p);

          std::ostringstream oss;
          oss << "p**" << boost::lexical_cast<std::string>( OrderChamp  ) ;
          funcs.push_back(oss.str());
        */
        /*    oss.str("");
              std::vector<std::string> funcs2;
              oss << "p**" << boost::lexical_cast<std::string>( Order ) ;
              funcs2.push_back(oss.str());
        */
        //test_app->addOutput( Output(_name="norm_L2",_latex="\\left\\| . \\right\\|_{L^2}",_dependencies=depend,_funcs=funcs) );
     
        //test_app->addOutputValue( 3.1415926 );
        //test_app->postProcessing();

        typedef Mesh<Simplex<Dim,OrderGeo,Dim> > mesh_type;
        typedef boost::shared_ptr<  mesh_type > mesh_ptrtype;

        typedef Mesh<Simplex<Dim,1,Dim> > mesh_1_type;
        typedef boost::shared_ptr< mesh_1_type > mesh_1_ptrtype;

        typedef bases<Lagrange<1,Scalar,PointSetFekete> > basis_1_type;
        typedef bases<Lagrange<2,Scalar,PointSetFekete> > basis_2_type;
        typedef bases<Lagrange<3,Scalar,PointSetFekete> > basis_3_type;
        typedef bases<Lagrange<4,Scalar,PointSetFekete> > basis_4_type;
        typedef bases<Lagrange<5,Scalar,PointSetFekete> > basis_5_type;
        typedef bases<Lagrange<6,Scalar,PointSetFekete> > basis_6_type;
        typedef bases<Lagrange<7,Scalar,PointSetFekete> > basis_7_type;
        typedef bases<Lagrange<8,Scalar,PointSetFekete> > basis_8_type;
        typedef bases<Lagrange<9,Scalar,PointSetFekete> > basis_9_type;

        typedef FunctionSpace<mesh_type, basis_1_type> space_1_type;
        typedef FunctionSpace<mesh_type, basis_2_type> space_2_type;
        typedef FunctionSpace<mesh_type, basis_3_type> space_3_type;
        typedef FunctionSpace<mesh_type, basis_4_type> space_4_type;
        typedef FunctionSpace<mesh_type, basis_5_type> space_5_type;
        typedef FunctionSpace<mesh_type, basis_6_type> space_6_type;
        typedef FunctionSpace<mesh_type, basis_7_type> space_7_type;
        typedef FunctionSpace<mesh_type, basis_8_type> space_8_type;
        typedef FunctionSpace<mesh_type, basis_9_type> space_9_type;

        typedef boost::shared_ptr<space_1_type> space_1_ptrtype;
        typedef boost::shared_ptr<space_2_type> space_2_ptrtype;
        typedef boost::shared_ptr<space_3_type> space_3_ptrtype;
        typedef boost::shared_ptr<space_4_type> space_4_ptrtype;
        typedef boost::shared_ptr<space_5_type> space_5_ptrtype;
        typedef boost::shared_ptr<space_6_type> space_6_ptrtype;
        typedef boost::shared_ptr<space_7_type> space_7_ptrtype;
        typedef boost::shared_ptr<space_8_type> space_8_ptrtype;
        typedef boost::shared_ptr<space_9_type> space_9_ptrtype;

        typedef typename space_1_type::element_type element_1_type;
        typedef typename space_2_type::element_type element_2_type;
        typedef typename space_3_type::element_type element_3_type;
        typedef typename space_4_type::element_type element_4_type;
        typedef typename space_5_type::element_type element_5_type;
        typedef typename space_6_type::element_type element_6_type;
        typedef typename space_7_type::element_type element_7_type;
        typedef typename space_8_type::element_type element_8_type;
        typedef typename space_9_type::element_type element_9_type;

        //-----------------------------------------------------------------------------------//
        double meshSize1=0.5;//test_app->vm()["hsize1"].template as<double>();
        //double meshSize1=0.7;//0.1;
        double meshSize2=1/64.;//256.;// 2/32.;//0.08;

        //-----------------------------------------------------------------------------------//
    
        //Mesh

        std::string __nameTest="convergence_p";
        //std::ostringstream __ostr;__ostr<<
        std::ostringstream __ostr_name1;
        __ostr_name1 << __nameTest << "_mesh_rect1_h" << meshSize1 <<"_";
        std::ostringstream __ostr_name2;
        __ostr_name2 << __nameTest << "mesh_rect2_h" << meshSize2 <<"_";
        std::ostringstream __ostr_name3;
        __ostr_name3 << __nameTest << "mesh_rectcirc1_h" << meshSize1 <<"_";
        std::ostringstream __ostr_name4;
        __ostr_name4 << __nameTest << "mesh_rectcirc2_h" << meshSize2 <<"_";
        //mesh_ptrtype mesh1 = createMesh_rect_p1<Dim,OrderGeo>( __ostr_name1.str(), meshSize1 );
        //mesh_ptrtype mesh2 = createMesh_rect_p2<Dim,OrderGeo>( __ostr_name2.str(), meshSize2 );
        //mesh_ptrtype mesh1 = createMesh_spec_p1<Dim,OrderGeo>( __ostr_name1.str(), meshSize1 );
        //mesh_ptrtype mesh2 = createMesh_spec_p2<Dim,OrderGeo>( __ostr_name2.str(), meshSize2 );

        mesh_ptrtype mesh1 = createMesh_rectcirc_p1<Dim,OrderGeo>(__ostr_name3.str(), meshSize1 );
        mesh_ptrtype mesh2 = createMesh_rectcirc_p2<Dim,OrderGeo>( __ostr_name4.str(), meshSize2/6. );
        //mesh_1_ptrtype mesh2 = createMesh_rectcirc_p2<Dim,1>( __ostr_name4.str(), meshSize2/6. );
    
        //-----------------------------------------------------------------------------------//
    
        space_1_ptrtype Xh1 = space_1_type::New( mesh1 );uint nb1dof=Xh1->nDof();
        space_2_ptrtype Xh2 = space_2_type::New( mesh1 );uint nb2dof=Xh2->nDof();
        space_3_ptrtype Xh3 = space_3_type::New( mesh1 );uint nb3dof=Xh3->nDof();
        space_4_ptrtype Xh4 = space_4_type::New( mesh1 );uint nb4dof=Xh4->nDof();
        space_5_ptrtype Xh5 = space_5_type::New( mesh1 );uint nb5dof=Xh5->nDof();
        space_6_ptrtype Xh6 = space_6_type::New( mesh1 );uint nb6dof=Xh6->nDof();
        space_7_ptrtype Xh7 = space_7_type::New( mesh1 );uint nb7dof=Xh7->nDof();
        space_8_ptrtype Xh8 = space_8_type::New( mesh1 );uint nb8dof=Xh8->nDof();
        space_9_ptrtype Xh9 = space_9_type::New( mesh1 );uint nb9dof=Xh9->nDof();
    
        element_1_type u1( Xh1, "u1" );
        element_2_type u2( Xh2, "u2" );
        element_3_type u3( Xh3, "u3" );
        element_4_type u4( Xh4, "u4" );
        element_5_type u5( Xh5, "u5" );
        element_6_type u6( Xh6, "u6" );
        element_7_type u7( Xh7, "u7" );
        element_8_type u8( Xh8, "u8" );
        element_9_type u9( Xh9, "u9" );

        //-----------------------------------------------------------------------------------//
        /*
          AUTO (h , cos(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
          AUTO (dxh , -0.5*M_PI*sin(0.5*M_PI*Px())*(cos(0.5*M_PI*Py())));
          AUTO (dyh , -0.5*M_PI*cos(0.5*M_PI*Px())*(sin(0.5*M_PI*Py())));
          AUTO (grad_h, vec(dxh,dyh));
        */
        AUTO( pi, M_PI );
        AUTO( h, cos(pi*Px()/2.)*cos(pi*Py()/2.) );
        AUTO( grad_h , vec( -(pi/2.)*sin(pi*Px()/2.)*cos(pi*Py()/2.) , -(pi/2.)*cos(pi*Px()/2.)*sin(pi*Py()/2.) ) );
    
        projectionH1<element_1_type>(u1,test_app);
        projectionH1<element_2_type>(u2,test_app);
        projectionH1<element_3_type>(u3,test_app);
        projectionH1<element_4_type>(u4,test_app);
        projectionH1<element_5_type>(u5,test_app);
        projectionH1<element_6_type>(u6,test_app);
        projectionH1<element_7_type>(u7,test_app);
        projectionH1<element_8_type>(u8,test_app);
        projectionH1<element_9_type>(u9,test_app);

        /*
          u1 = project( Xh1, elements(mesh1), h );
          u2 = project( Xh2, elements(mesh1), h );
          u3 = project( Xh3, elements(mesh1), h );
          u4 = project( Xh4, elements(mesh1), h );
          u5 = project( Xh5, elements(mesh1), h );
          u6 = project( Xh6, elements(mesh1), h );
          u7 = project( Xh7, elements(mesh1), h );
        */

        //-----------------------------------------------------------------------------------//

        std::ostringstream __ostr;
        __ostr<<"graph_boundary_conv_p_"<<OrderGeo<<".data";
        std::ofstream out;
        out.open(__ostr.str().c_str(), std::ios::out);

        //double  __aera = integrate( markedfaces(mesh2, mesh1->markerName("Interface")),cst(1.0) ).evaluate()(0,0);
        double  __aera = integrate( markedfaces(mesh2, mesh1->markerName("Interface")), idv(u1)*idv(u1) ).evaluate()(0,0);

        //std::vector<boost::tuple<double,double> > __err(9);

        integration<element_1_type,mesh_ptrtype,OrderGeo>(u1,mesh2,out);
        integration<element_2_type,mesh_ptrtype,OrderGeo>(u2,mesh2,out);
        integration<element_3_type,mesh_ptrtype,OrderGeo>(u3,mesh2,out);
        integration<element_4_type,mesh_ptrtype,OrderGeo>(u4,mesh2,out);
        integration<element_5_type,mesh_ptrtype,OrderGeo>(u5,mesh2,out);
        integration<element_6_type,mesh_ptrtype,OrderGeo>(u6,mesh2,out);
        integration<element_7_type,mesh_ptrtype,OrderGeo>(u7,mesh2,out);
        integration<element_8_type,mesh_ptrtype,OrderGeo>(u8,mesh2,out);
        integration<element_9_type,mesh_ptrtype,OrderGeo>(u9,mesh2,out);


    } //end run_convergence_p





}//end namespace test_interp_twomesh



    /*_________________________________________________*
     *_________________________________________________*
     *_________________________________________________*
     * MAIN
     *_________________________________________________*
     *_________________________________________________*
     *_________________________________________________*/



    //BOOST_AUTO_TEST_SUITE( interp_twomesh_testsuite )

    //BOOST_AUTO_TEST_CASE( interp_twomesh )
    int main(int argc,char** argv)
    {

        using namespace test_interp_twomesh;

        Application_ptrtype test_app(new Application_type(argc,//boost::unit_test::framework::master_test_suite().argc,
                                                          argv,//boost::unit_test::framework::master_test_suite().argv,
                                                          makeAbout(),
                                                          makeOptions()
                                                          ));

        test_app->changeRepository( boost::format( "/testsuite/feeldiscr/%1%/" )
                                    % test_app->about().appName()
                                    );




        //run_test_geomap<2,4,1>(test_app);
        //run_test_element<2,4,1>(test_app);
        //run_test_element<2,4,2>(test_app);
        //run_test_element<2,4,3>(test_app);

        /*
          Log() << "\n[main] ----------TEST_BOUNDARY_START----------\n";
          Log() << "[main] ----------------<2,4,1>---------------\n";
          run_test_boundary<2,4,1>(test_app);
          Log() << "[main] ----------------<2,4,2>---------------\n";
          run_test_boundary<2,4,2>(test_app);
          Log() << "[main] ----------------<2,4,3>---------------\n";
          run_test_boundary<2,4,3>(test_app);
          Log() << "[main] ----------------<2,4,4>---------------\n";
          run_test_boundary<2,4,4>(test_app);
          Log() << "[main] ----------------<2,4,5>---------------\n";
          run_test_boundary<2,4,5>(test_app);
          Log() << "[main] ----------TEST_BOUNDARY_FINISH----------\n\n";
        */



        test_app->changeRepository( boost::format( "/Test_Conv_Interp/%1%/" )
                                    % test_app->about().appName()
                                    );

        //run_test_convergence_h<2,1,1>(test_app);
        //run_test_convergence_p<2,4,1>(test_app);

        run_test_convergence_h<2,4,2>(test_app);
        //run_test_convergence_p<2,4,2>(test_app);



        //run_test_boundary2<2,4>(test_app);

        Log() << "[main] ----------------<2,6,2>---------------\n";
        //run_test_geomap<2,4,2>(test_app);
        Log() << "[main] ----------------<2,6,3>---------------\n";
        //run_test_geomap<2,2,3>(test_app);
        /*Log() << "[main] ----------------<2,6,4>---------------\n";
          run_test_geomap<2,6,4>(test_app);
          Log() << "[main] ----------------<2,6,5>---------------\n";
          run_test_geomap<2,6,5>(test_app);*/
        Log() << "[main] ----------TEST_GEOMAP_FINISH----------\n\n";
        Log() << "[main] ----------TEST_INTERP_START-----------\n";
        /*    Log() << "[main] ----------------<2,7,1>---------------\n";
              run_test_interp<2,7,1>(test_app);
              Log() << "[main] ----------------<2,8,2>---------------\n";
              run_test_interp<2,8,2>(test_app);
              Log() << "[main] ----------------<2,9,3>---------------\n";
              run_test_interp<2,9,3>(test_app);
              Log() << "[main] ----------------<2,10,4>---------------\n";
              run_test_interp<2,10,4>(test_app);
              Log() << "[main] ----------------<2,11,5>---------------\n";
              run_test_interp<2,11,5>(test_app);
              Log() << "[main] ----------TEST_INTERP_FINISH----------\n";
              Log() << "[main] ----------TEST_EXPORT_START-----------\n";
              Log() << "[main] ----------------<2,7,1>---------------\n";
              run_test_export<2,7,1>(test_app);
              Log() << "[main] ----------------<2,8,2>---------------\n";
              //run_test_export<2,8,2>(test_app);
              Log() << "[main] ----------------<2,9,3>---------------\n";
              //run_test_export<2,9,3>(test_app);
              Log() << "[main] ----------TEST_EXPORT_FINISH----------\n";*/

        return 0;
    }
//BOOST_AUTO_TEST_SUITE_END()
